---
title: "Figure 3"
output: html_document
---

```{r setup, include=FALSE}
setwd("/Users/shanhe/Desktop")
```

```{r}
r <-seq(0.02,0.05,0.01)
rho <-seq(0.01,0.04,0.01)
delta <-seq(0.7,1.0,0.1)
t<-seq (1,60,1)
# create all combinations 
all_cbs<- expand.grid(r, rho, delta,t) 
```

```{r}
#rename variables
colnames(all_cbs)
names(all_cbs)[names(all_cbs) == "Var1"] <- "r"
names(all_cbs)[names(all_cbs) == "Var2"] <- "rho"
names(all_cbs)[names(all_cbs) == "Var3"] <- "delta"
names(all_cbs)[names(all_cbs) == "Var4"] <- "t"
```
```{r}
all_cbs$P=2.981+0.328*(all_cbs$t)
all_cbs$P_star=2.981+0.29+0.328*(all_cbs$t)
all_cbs$D=0
all_cbs$D_star=0
```

```{r}
all_cbs$tau= (all_cbs$rho/all_cbs$r)*((1-exp(-all_cbs$r*all_cbs$delta*all_cbs$t))/(1-exp(-all_cbs$rho*all_cbs$t)))*(all_cbs$P_star+exp(-1*all_cbs$rho*all_cbs$t)*all_cbs$D_star)-(all_cbs$P+exp(-1*all_cbs$r*all_cbs$delta*all_cbs$t)*all_cbs$D_star)
```

```{r}
#get data s.t r>rho
sdata <- subset(all_cbs, all_cbs$r>=all_cbs$rho, select=c(r,rho, delta, t, P, P_star, D, D_star, tau))
summary(sdata$tau)
plot(sdata$t, sdata$tau)
write.table(sdata, file="figure_3.csv",sep=",",row.names=F)
```